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Wavelet-based Estimator for the Hurst 
Parameters of Fraetional Brownian Sheet 


Liang Wu and Yiming Ding 


Abstract 

It is proposed a class of statistical estimators H = {Hi ,... ,Hd) for the Hurst parameters H = 
(Hi,..., Hd) of fractional Brownian field via multi-dimensional wavelet analysis and least squares, 
which are asymptotically normal. These estimators can be used to detect self-similarity and long-range 
dependence in multi-dimensional signals, which is important in texture classification and improvement 
of diffusion tensor imaging (DTI) of nuclear magnetic resonance (NMR). Some fractional Brownian 
sheets will be simulated and the simulated data are used to validate these estimators. We find that when 
Hi > 1/2, the estimators are efficient, and when Hi < 1/2, there are some bias. 

Index Terms 

Detection of long-range dependence. Self-similarity, Hurst parameters. Wavelet analysis. Fractional 
Brownian sheet 


I. Introduction 

In numerous fields such as hydrology, biology, telecommunication and economics, the data available 
for analysis usually have scaling behavior (or long-range dependence and self-similarity) that need to be 
detected. The key point for detecting scaling behavior is the estimation of the Hurst parameters H, which 
are used in characterizing self-similarity and long-range dependence. The literature on the topic is vast 
(see e.g. Hl-S). Most of the published research concerns the scaling behavior in ID data. However, 
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it is also important to detect the scaling behavior in multi-dimensional signals, which is important in 
texture classification and improvement of diffusion tensor imaging (DTI) of nuclear magnetic resonance 
(NMR). Moreover, many multidimensional data from various scientific areas also have anisofropic nafure 
in the sense that they have different geometric and probabilistic along different directions. The study 
of long-range dependence and self-similarity in multidimensional case is usually based on the model of 
anisotropy fractional Gaussian fields such as fractional Brownian sheet (fBs) ifTOl . ifTTI . operator scaling 
Gaussian random field (OSGRF) ifT^ and extended fractional Brownian field (EFBF) ifT^ . 

We focus on fBs, which plays an important role in modeling anisotropy random fields wifh self- 
similarity and long-range dependence (see e.g. ifTTl . O, ifTSl l since it was first introduced by Kamont 
ifTOl . Fractional Brownian sheet has arisen in the study of texture analysis for classification ifT^ . ifTTl . 
Besides research on the texture analysis, fractional Brownian sheets can also be used to drive stochastic 
partial differential equations ifTSl . Moreover, Diffusion-tensor images (DTI) of nuclear magnetic resonance 
(NMR) in brain is now based on three-dimensional Brownian motion |[T^ Hi = H 2 = H^. Compared 
with this model, fractional Brownian sheet may be a better model because it can describe the long-range 
dependence, which may exist in brain. 

The main purpose of this paper is to estimate the Hurst parameters of fractional Brownian sheets 
via wavelet analysis. Wavelet analysis has been an efficient tool for the estimation of Hurst parameter, 
but most of the published research on the topic focus on the ID case. Wavelet-based estimators for the 
parameter of single-parameter scaling processes (self-similarity and long-range dependence) have been 
proposed and studied by Abry et al. (see Hl-lHl). Such estimators are based on the wavelet expansion of 
the scaling process and a regression on the log-variance of wavelet coefficient. It was usually assumed that 
wavelet coefficients are uncorrelated within and across octaves. Bardet ll^ and Morales ll^ removed this 
assumption and proved the asymptotic normality of these estimators in the case of fractional Brownian 
motion (fBm), proposed a two-step estimator, which has the lowest variance in all of the least square 
estimators. Compared with other estimators, such as the R/S method, the periodogram, the maximum- 
likelihood estimation, and so on, the wavelet-based estimator performs better in both the statistical and 
computational senses, and is superior in robustness (see S-llhl and the references therein). It not only 
has small bias and low variance but also leads to a simple, low-cost, scalable algorithm. Besides, the 
wavelet-based method can also eliminate some trends (linear trends, polynomial trend, or more) by the 
property of its vanishing moments 0, which makes the estimator robust to some nonstationarities. 

The paper is organized as follows. In Section ini some properties of Gaussian random variables and 
some limiting theorems used in this paper are introduced. In Section |III1 we obtain some properties of 
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wavelet coefficient of fBs (see Proposition [T]). Based on these properties, A class of estimators H for 
the Hurst parameter of fBs are proposed (see Section ITTT-BIi . We also prove that these estimators H are 
asymptotically normal following the approach in fl\\ (see Theorem |4] and Section ITTT-Ch . Using the two- 
step procedure that proposed hy Bardet ll22l . we realize the generalized least squares (GLS) estimator, 
which has the lowest variance of H (see Section ITTT-Dh . Finally, this two-step estimator Hog is applied 
to the synthetic fBs generated hy the method of circulant embedding in the two-dimensional case and 
three-dimensional case (see Section HVl) . 

II. Preliminaries 

Definition 1. Given a vector H = {Hi,..., H^) E (0,1)'^, a fractional Brownian sheet (fBs) = 
{B^{t),t G with Hurst parameter H is a real-valued mean-zero Gaussian random fields with 
covariance function given by 

E{B^{s)B^{t)) = ^ - I'Si - G M+. (1) 

i=l 

It follows that B^ is an anisotropic Gaussian random field. 

Note that if d = 1, then B^ is a fractional Brownian motion (fBm) with Hurst parameter Hi G(0,1), 
which as a moving-average Gaussian process is introduced hy Mandelbrot and Van Ness E^ . If d > 1 
and Hi = ■ ■ ■ = H^ = 1/2, then B^ is the Brownian sheet. 

Moreover, B^ has a continuous version for all H = (Hi ,..., Hfj G (0,1)'^ and is self-similar in the 
sense that for any diagonal matrix A = diag(ai, 02 , ■ ■ ■, a^), where Oj > 0 for all 1 < i < d, the random 
field {B^{Af),t G has the same probability law as {a^B^{t),t G M'^}, where 
(see e.g. my B^ also has stationary increments with respect to each variable (see e.g. EH). 

First we will use some properties of Gaussian random variables and some limiting theorems. 

Lemma 1. {X,Y) are two-dimensional mean-zero normal variable, i.e., {X,Y) N{0,al-,0,al-,p), 

then 

cov{X^,Y'^) =2cov^{X,Y). 

Proof: Assume cov{X,Y) = v. Let Z = aX — Y s.t. cov{Z,Y) = 0, So we have a = 

= o^EX^ Ey2 _ 2aEXy = a^a\ Pal- 2av. 

EX2y2 ^ J_]Ez2y2 + ^EZy3 + ^rEy^ 

H a H 

= ^rEz^y^ -EZEy^ + ^rEy"^. 

a2 a a2 
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It’s easy to check = 0. Put a and EZ^ into the above formula, we have EX^Z^ = (rfirl + 2u^. 

Then cov(X^,y^) = 2cov^(X, Z). ■ 

Definition 2. Let f ^ R. X = .. ,XW)' is a Gaussian vector, and f{X) has finite second 

moment, we define the Hermite rank of f with respect to X as 

rank(/) := inf{T : 3 ploynomial P of degree r with E[(/(X) — Ef{X))P{X^^\ ... / 0}. 

The following two lemmas have been proved in ll^ (see Lemma 4 and Lemma 5 of ||2T1). 

Lemma 2. If X ^ N{Q,a‘^), then for q > —1/2, the rank of fq{X) = \X\‘t is 2. 

Lemma 3. Let hfix) = fq{x) = \x\^,q > —\I2 for every 1 < i < d Let ti,... fid be d real numbers, 
and X = a zero-mean Gaussian vector. Then 

h{x) = yf, tMxy 

is a function g : —)• R o/ Hermite rank 2. 


The following three limiting theorems have been proved in |[26l . lIlTl and ll28l (see Theorem 9.0 of 
l|26l, Theorem 4 of ll27l and Theorem 3.3A of ll^Sl t. 


Theorem 1 (Cramer-Wold theorem). Let Xn = {X^^^,..., Xlf^) and X = {X^'^\...,Xy be random 
vectors. As n ^ +oo, Xn X if and only if: 




for each {ti,... fid) G R'^. 


Theorem 2 (Arcones, 1994). Let {Xj}^^ be a stationary mean-zero Gaussian sequence of R'^-valued 
vectors. Set Xj = {Xj^\ ..., Let f be a function on R'^ with rank r, 1 < r < oo. We define 

for /c G Z, where m is any large enough number such that k > 1. Suppose that 

y <oo, ( 2 ) 

for each 1 < iifi 2 ^ d. Then 
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where 


Theorem 3. Let Xn = {Xn'^ 


= Var/(Xi) + 2 Cov{f{X,), f{X,+k))- 

^-^k=l 

, xif^). Suppose that as n ^ oo, 


( 4 ) 


Xn p d 


>iV(0,S), 


where S is covariance matrix and bn —)• oo. Let g{x) = {gi{x),..., gm{x))■, x = (xi,..., x^), be 
a vector-valued function for which each component function gi{x) is real-valued and has a nonzero 
differential gi{p;t),t = (ii,..., td), at x = p.Put 

' dgi . 


D = 


dxr=^. 


mxd 


Then 


g{Xn) - g{p) d 


n{o,dj:d') 


III. Wavelet-based estimator 
A. Wavelet coefficients of fractional Brownian sheet 

Let 'f{f) G L^(]R) be an orthonormal wavelet with compact support and have > 1 vanishing 
moments. A wavelet is said to have a number > 1 of vanishing moments if 


/ 


t^'f{t)dt = 0, fc = 0,1, 2,... , A^ — 1, and / t^Tp{t)dt / 0. 


The dilations and translations of '^{t) 

-k), j,k 


(5) 


( 6 ) 


construct the wavelet orthonormal basis for L^(R). The factors 2^ and j are called the scale and octave 
respectively. The factor k is called translation. The wavelet orthonormal basis for can be obtained 

by simply taking the tensor product functions generated by d one-dimensional bases: |[29l . Il30l 


■ ■ ■ fd) — (^l)V'i2A2 (^ 2 ) • • • '4’jd,kd{td)- 


(V) 


For convenience, let J = (ii,... ,jd) £ K = (fei, ..., kd) G t = {ti,... ,td) G Mf. Then the 
wavelet coefficients of fractional Brownian sheet in are: 

dBH{J,K):= [ ■■■ [ 'fj^^kAti)'^j2M{t^)---'^3ci,kd{td)B^{t)dti---dtd. (8) 

dR Jm. 

Before introducing our estimator for , we present some properties of the wavelet coefficient of B^ . 
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Proposition 1. Let ijjlt) G L^(M) be an orthonormal with compact support and have N > 1 vanishing 
moments. The wavelet coefficients of fractional Brownian sheet defined by Eq.^ have the following 
properties: 

(i) KdsH {J, K) = 0, and dBH{J,K) is Gaussian, for any J,K £ 

(ii) For fixed J = ■ ■ ■ ,jd) £ tind \/Ki ,..., Kn G we have 

{dBH{J, Kl), . . . , dBH{J, Kn)) = Ki), . . . , dBH{0, Kn)), (9) 

where 0 denotes (0,0, • • • ,0) G 

(in) For fixed J G Z'^, and Vi^i,..., Kn, L G Z'^, we have 

{dBH{J,K^ +L),.. .,dBH{J,Kn + L)) ={dBH{J,K^), . . .,dBH{J,Kn)), (10) 

(iv) For fixed J = (ji,... ,jd) G J' = (j'l, • • • ,jd) G Z'^, K = (fci,..., G Z"^, and K' = 

(k[,... ,k'f) G Z'^, when min \2^^ki — +oo 

l<i<d 

d 

EdBHiJ,K)dBH{J',K') « Yl\2^^ki - 2^'k'f^'-^^. (11) 

i=l 

In the above, = means equality in distribution. 

Proof: (i) 

mBH{J,K)= [ EB^{t)i:j^^k^iti)'fj,,k2{t2)---'ipu,kMdt, 

EB^{t) is a constant, and f 'fj^k{t)dt = 0, so EdBH{J,K) = 0. The Gaussianity of dBH{J,K) follows 
from the fact that the wavelet coefficient of Gaussian process is Gaussian (see |[3ll ). 

(ii) This property can he obtained hy the self-similarity of B^. For MK G 

dBH{J,K) = / V'(2-^Hi - ki) ■ • • V’(2"^"fd - kd)B^{t)dt 

Jr'^ 

= 22^-1^* / 'f{t'i-ki)---'f{t'd-kd)B^{2^H[,--- ,2^H'f}dt' 

= dBH{^,K), 

where t' = (t'l,..., t'f). 

Similarly, one can check that, for G M, 1 < i < n, 

V" = 6idBn{Q,Ki). 

z—/j=l ^ —'1=1 
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(iii) Let d{J, ki,..., kd) = d{J, K), L = (^i,..., Id). Without loss of generality, we can take J = 0. 
First we check the formula helow, 

ds^iO, ki + h,... ,kd + ld) = dB^iO, ki,k2 + h, ■ ■ ■ ,kd + Id)- ( 12 ) 

In fact, 


dBH {0,ki + li,... ,kd + Id) 

= [ ipiti -ki-h)--- 'ijjitd - kd- ld)B^{ti,td)dti ■■■dtd 


= / - ki)B^{t'i td)dt'i ■ IT ip{ti - ki- li)dt2 ■ ■ ■ dtd- 

jR-i-i Jr ^^2 

Since — ki)dt'i = 0 and B^ has stationary increments with respect to each variable, we have 

[ 'il;{t[-ki)B^{t[+li,t2...,td)dt[ = [ 'il;{t[-ki)B^{t[,t2...,td)dt[. 


The Eq.lfT^ is proved. 

Applying Eq. dT^ repeatedly, 

dBH (0, K + L) = dBH (0, fci, /c2 + I 2 , ■ ■ ■, kd + Id) — dBH (0, ki, k 2 , k^ + I 3 ,..., kd Id) 
= ... =dBH (0, /ci,..., kd) = dBH{0, K). 


Similarly, one can see that, for G M, 1 < i < n, 

E TT. (] » ^ 7T. 

, + = eidBH{0,Ki). 

1 = 1 /2 = 1 

(iv) The covariance function of B^ is given hy Eq.©. And 


mBHiJ,K)dBH{f,K' 


KB^ {t)B^ { 3 )^ 111 ^ji,kA^i)'il^fi,k'{si)dsdt. 


2=1 


So we have: 


EdBH{J,K)dBH{J',K') 

= n / / - |si - ti\^^')dsidti 

Jr Jr ^ 

d 

= EdBHi {ji, ki)dBHi {j'i, k'i), (13) 

2=1 

where B^' is fBm with Hurst exponent Hi. 
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It has been shown that the correlations between wavelet coefficients of fBm satisfy this asymptotic 
equation ||6l, ||32l, 1331 

EdBH{j,k)dBH{j',k') « \2^k — 2^ k'\'^^~‘^^, \2^k — 2^ k'\ —>• oo, 

where j,/, k, k' G Z, and is fBm with Hurst exponent H, and the wavelet used has compact support. 
The Eq. dmi is proved. ■ 

Remark 1. In view of Eq.^n\. to avoid long-range dependence for dBH{J,K), i.e., to ensure that 
^\dB^{Ji K)dB^ ifli 0)1 < oo, one can choose 

N> max i7j + 1/2, i.e., max (2Tfj — 2A^) < —1. (14) 

l<i<d l<i<d 

Under this condition, the correlation of dBH{J,K) tends rapidly to 0 at large lags. 

B. Estimation of the Hurst parameter and its asymptotic behavior 

In this subsection, we estimate the Hurst parameter of fractional Brownian sheet in the d-dimensional 
finite interval. The following lemma can be obtained by direct calculation. 

Lemma 4. G ~ {Ti^---)Td) £ is a fractional Brownian sheet in the 

d-dimensional finite interval, ff) is an orthogonal wavelet with compact support [— M, M], M > 0 and 
has N > 1 vanishing moments. We define: 

dB^,TiJ,K) ■■= [ ■■■[ fj^,k^itl)'fj2,kAt2)■■■'^l:j,,,kAtd)B^ii)dtl■■■dtd■ (15) 

Jo Jo 

For ease of writing, let M <\. So when 1 <ki < 2~^*Ti — M, z = 1, ..., d, 

dB»,T{J,K)=dBn{J,K), (16) 

where dBH{J,K) is defined by Eq.^. 

Definition 3. Under the conditions of Lemma ^ we define that the available wavelet coefficients for the 
finite interval fractional Brownian sheet are the wavelet coefficients Jbh K) that satisfy Eq.l\16\). 
Let rij be the number of available wavelet coefficients at octave J G Z'^. 

nj = n rii, Hi := [2“-^'*Ti - MJ . 

If there exists zq G {1,... , d} s.t. Ti^ ^ +oo, then nj +oo. 
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Under the conditions of Lemma |4l according to Eq.@, Eq. lfTOl) and EemmalU the wavelet coefficients 
of (t) satisfy the equation, 


EdlH^TiJ,K) =e4k(J,K) = (7^^2Ji(2/ri+i)+j2(2H.+i)+...+i,(2/r,+i)^ 


(17) 


where the Hurst parameter H = {Hi,..., H^) G (0,1)'^, Cb^ = Ed^H(0,l) only depends on , 

J = (ji,- • • ,jd),K = {ki,.. .,kd)G 

Take the 2-logarithm of both sides: 

log2 Ed^H J, K) = yi(277i + 1) + • • • + jd{2Hfi + 1) + log2 Cb^ . (18) 

So the estimation of H can he realized hy a linear regression in the left part versus J diagram. The 
Ed^K rp{J,K) can he estimated hy 

S(J) 

According to Eemma |4l 

s(j)=i/n,E:;,-Ei;.,4-«/f). 

In more detail, the H is estimated hy the regression 

Ji 1 ^ 


(19) 


/ log 2 S(Ji) \ / 


/ 2idi + 1 \ 


y iog2 5’(Jm) / 


V 


Jm. 1 


/ 


2Hci + 1 

y iog2C'BH y 


/ e(Ji) \ 


y £{Jm) J 


where Ji = {ji,i, ■ ■ ■ ,ji,d) G G m}, and Ji < • • • < Jm, { Ji < Jk in the sense that 

ji,i < jk,i, for i G {1,... , d} and the equality can’t hold for all i ). 

Eet 

/ 2idi + l \ 


L := 


( log2 5(Ji) ^ 


y log 2 5 ’(Jm) / 


,A:= 


/ 


V 


Ji 1 


Jm. 1 


\ 


/ 


,a := 


2Hd + 1 

y iog2 Cbh y 


,e:= 


/ e(Ji) \ 


y ^{Jm) J 


suppose that m > d + 1, rank(A) = d + 1. 

Then least squares estimation gives the estimator: 

d = (A'E-^A)-^AE-^L, 

H = a/2- 1/2, 


( 20 ) 
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where H = (^i, • • • , H^, log 2 Cb ^/2 — 1 / 2 )'. Hi is the wavelet-hased estimator of Hi, i G { 1 ,... , d}, 
S is a m X m full-rank matrix. H is selected to minimize 

\\L - Aa||| = (L - AayT,-^{L - Aa). 

If S = Imxm, H is the ordinary least squares (OLS) estimator. If S = VarL = Vare, H is the 
generalized least squares (GLS) estimator. 

The following central limit theorem indicates the asymptotic normality of the estimators H. 

Theorem 4. Under the conditions of Lemma^ H is defined by Eg. d20l). N > max Hi + 1/4. For any 

l<i<d 

c G { 1 ,..., d}, Tc ^ +00, so ni^c +oo and nj^ +oo, / G { 1 ,... , m} (nj^ = nf=i same 

as that in Definition |21), we have 

nUUH-U)ANiO,J:^), (21) 

where U = (Tf', log 2 C'bh/2 — 1/2)' and 

( 22 ) 

with Tjb = D'TjsD, S 5 = {crfjfjmxmy ^“ih defined by and 

D = log 2 e • diag( [min (e4h {Jm,K))~^). 

Remark 2. Tiijnm,,c i^ the asymptotic covariance of L. For GLS estimator, S is the covariance matrix 
of L. IfTc^ +00, S —^ Si/ nm,c> then 

S^ = ^{A'J:yA)-\ (23) 

Remark 3. According to the equation between Tc and nm,c (Definition 13), Theorem |4] indicates that 
the convergence rate of H is 0{l/^/Tf) as Tc —)• +00 for any c G {1,... , d}. Similarly, we can also 
obtain that the convergence rate of H is 0{1/\J ]/[^=i '^i) +00 for all i G {1,... , d}. That is 

to say, the estimation accuracy of H can be improved by increasing the square root of sample volume 

\/Y\Un 

C. Proof of Theorem |4] 

The following limit theorem of S{J) can he obtained hy using Theorem |2] 
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Lemma 5. Under the conditions of Lemma^ S{J) is defined by EqAl%. N > max Hi + \/A. For any 

l<i<d 


c G {1,... ,d}, Tc ^ + 00 , so Tic —)• +00 and nj —>• +oo, we have 

nll^[S{J)-u{J)]ANf),a\j)), 
where nj is the same as that in Definition\^ u(J) = {J, K), 


(24) 


a^iJ) = 


{njjncY 


Var( dl„{J,K)) 


kc=l^kiG{l,...ni},i^c 


+ 


'' ■^1 k=l kc=l,kie{l,--ni},i^c 


dl^{J,K)). (25) 


kc=l-\-k,ki£{l,...ni},i^c 

Proof: We construct a sequence of random vectors Xy from the wavelet coefficients dsH {J, K), 


Xy = {Xyil),...,Xy{nj/nc)y 


I fce=D,fci(nj/nc) X1) ^ F: V F: tic- 


When Ur 


+ 00 , using Proposition [TJ is a stationary mean-zero Gaussian sequence. Define 

1 


h{Xy) = 


nj/n, 


E nj/ric 2 

. , \Xv{i)V. 

2=1 


According fo LemmaO fhe Hermife rank of h{Xy) is 2. So hy Proposition [TJ when N > max Hi + l/A, 

l<i<d 

for each 1 < ii,f 2 < nj/uc. 


k(*i’* 2 )(/c)p < oo, 

where r(*i’*=)(A:) = Wf<:rn{k)Xm+k{i 2 )]- 

Hence, hy Theorem |2l fhe lemma is frue. ■ 

Now consider fhe asymptotic behavior of fhe full vector S = (5(Ji),..., S{Jm))'■ In order fo apply 
Theorem [U (Cramer-Wold fheorem), we need to study the asymptotic behavior of aS, a = (ai,... ,am) 
is a fixed buf arbifrary element of M”*. 


Lemma 6. Under the conditions of Lemma ^ let a = (oi,..., am) be a fixed but arbitrary element of 
M™, S = (S'( Ji),..., S{Jm))'- X > max Tfj-|-l/4. For any c € {1,... , d}, Tc —>■ -|-oo, nj, = n^,=i 

l<i<d ’ 

u{Ji) = (J/, K), f € {1,..., m}, we have 

nllliaS - J]]™ ^ aiu{Ji)] 4 iV(0, alg), (26) 
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Ks = Var 


2^Cov 

k=l 


\ 


E 

1=1 


ai- 


ni,c 


dBH{Ji,K) 


kie{l,...,ni^i},i^c 


+ 


\ 


E 


ai 


ni,c 


1 i. 1 


1=1 

fciG{l,...,ni,i},i^c 


ni. 


yy d%H{Ji,K) 


kc&{uik+l,...,ui{k+l)}, 


/ 


(27) 


aiiJ K = {ki,..., kd), ui = 2'^™''= m}. 


Proof: We construct a sequence of random vectors ^ from the wavelet coefficients dBH{Ji, K), 
K = {ki,...,kd). 


^l,s (-^l,s (1)) • • • ) ^l,s Ji/^l,c)) • ('7l) ) |fc^=s,fciG{l,,..ni_i},l^c] 1X (n jj/ni_c) ’ 1 E S E 

Because Ji < • • • < Jm, according to Definition [3l ni^c > • • • > Bm,c- One can check ui « Without 
loss of generality, we suppose that ni^c = 'u-i'^m^c, I £ {Ij ■ ■ ■ j "i}- 
Let 

= (^l,ui(t;—1)+1) ■ ■ ■ )^l,uit!) ■ ■ ■ ) 1)+1 J ' ' ' ^ ^m,u^v') j 1 E H E nm^c- 


When Tc —>• +oo, nm,c +oo, using Proposition [T] is a stationary mean-zero Gaussian 

sequence. Define 


h{x,) = 




By Lemma [3 the Hermite rank of h{X^) is 2. 
We rewrite X,, as 




Ul{v-l)+j 


E ttl 

So hy Proposition m when N > max L7j-|-l/4, it’s easy to check, for each 1 E * 1 , *2 E X]i=i ^i^Ji/^i,c> 

l<i<d ’ 

Y 7 < 00, 

^/c =—00 

where = E[Xjn{h)Xm+k{i 2 )]- 

Hence, hy Theorem |2j Lemma is proved. ■ 

So hy Theorem [T] (Cramer-Wold theorem), the asymptotic behavior of the vector S can he obtained. 
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Lemma 7, Under the conditions of Lemma^ S = {S{Ji),, S'( Jm))'. N > max Hi + 1/4. For any 

l<i<d 

c G {1,... ,d}, Tc ^ + 00 , we have 


nUl[S-u]ANiO,J:s), 


(28) 


where 


u — (Edg// (Ji, K ),..., Ed|jff (Jm, 7^)), 


(29) 


and 


with 


^2 _ 
^i.h — 


E 5 — {(^l^h)'mxmi 


uinjjni^c Uhnj^/rih^, 


■[Cov 


(30) 




5] dl„iJi,K), dlH{A,K) 




kc&{l,...,Uh}, 

kie{l,...,nh,i},ij^c 


+ ^ Cov 


k=l 


\ 




dlH(.Jh,K) 


kce{i,...,ui}, 

\ kiG{l,...,ni^i},i^c 

/ 


kce{uhk+l,...,Uh(k+l)}, 


+ ^ Cov 


k=l 


\ 


Y, dlu{JuK), Y d%n{Jh,K) 

kc;e{uik+l,...,ui{k+l)}, fccG{l,...,tth}, 

kie{l,...,nh,i},ij^c / 


ft/c t \ Lij n/-fl,.. 

\ki£{l,...,niH,i¥^c 

1 < l,h < m, and ui = nj^ = ^ € { 1 ,... , m}. 


(31) 


can be calculated by < 7 ^^ in Lemma 
Note that u/; is equal to a'^{Ji)/ui in Lemma |5] In fact, 


/ 


a?,= 




[Var 




E dl„{JuK) 


kce{i,...,ui}, 




/ 


+ 2 Yj Cov 


k=l 


Y dlu{JuK), Y 


\ 


dlniJuK) 


ki£{l,...,niH,i¥^c 


kc€{uik+l,...,ui{k+l)}, 


/ 




E Var 


E d%H{Ji,K) 
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j=i k=i 


d%H{JuK), ^ d^niJiiK) \]. 

I fcc=j,fciG{l,...,ni,i},i^c fce=j-|-fc,fciG{l,...,ni,i},i^c 


Since the sequence of random vectors in Lemma [5] is stationary, 
1 


= 


ui{njjni^c 


OO 


+ 2 ^ Cov 


k=l 
_ 2 / 


r[Var( d%n{Ji,K) 


.kc=l,ki£{l,...,ni^i},i^c kc=l+k,ki€^{l,...,ni^i},i^c 


= o- {Ji)/ui. 

By Theorem |3 one gets the asymptotic behavior of the logarithm of S. 


Lemma 8. Under the conditions of Lemma H] L = (log2 <S'( Ji),..., log2 S{Jm))'- > max Hi + 1/4. 

l<i<d 

For any c € {1,..., d}, Tc —)• +oo, we have 

n]ll[L-Lu]AN{Q,i:L), (32) 

where 

Lu = (log2lEd|H(Ji,iT),... ,log2lEd|H(Jm,-?^)), (33) 

and 

Si = D'l^s^D. (34) 

with D = log2 e • diag((Ed^ff (Ji, itT)) ^, • • •, (Ed^n {Jm,K)) ^), and S5 is defined by ddO]) and lU/l) . 
the same as that in Lemma 0 

The estimator H denoted hy Eq. (l20l) is linear combinations of the elements of L. So the asymptotic 
normality of H follows from the asymptotic normality of L by Theorem [3 


D. Two-Step estimator 

We adopt the GLS estimator, which has the minimum variance in these estimators defined by (l20l ) 
(see Gauss-Markov theorem |[34l l. But the covariance matrix VarL used by GLS estimator is a function 
of the unknown parameter H we want to estimate. So a two-step estimator may be feasible, namely, we 
first obtain the OLS estimator of H, and then use the OLS estimator to estimate the covariance matrix. 
So the GLS estimator can be obtained using the estimate of the covariance matrix. This estimator has 


Februaiy 4, 2015 


DRAFT 



15 


been used to the estimation of Hurst parameter of fBm and is proved to be asymptotically normal (see 

ma-ca). 


In this subsection, we apply two-step estimator to the case of fractional Brownian sheet. 

First we obtain the covariance matrix of L with known H. By the regularity on of logarithm 
function and Theorem |3 when nj^,nji is large, the {h,l) element of the covariance matrix VarL is 
given by 


Cov(log2 log2 S{Jl)) 


(log2e)"Cov(5(Jfe),g(J;)) 

mln{Jh,K)Edln{JuK) 


(log 2 e) 


E E Cov(4,(4, Kh)^ 4. {Jum 

2 Ki Kh _ 

nj^nj,E4 h {Jh, K)E(f^H {Ji,K) 


E 2Cov^((iB»(-^/i) Kh),dBH{Ji, Ki)) 


{log^e) 


2 Ki Kh 


nj^nj^V.d\H{Jh,K)^d\H{JhK) 


.(by Lemma [Hi 


(35) 


Consider Cov{dBH{Jh, Kh), ds^ {Ji, Ki)), by Eq.([l3]l, 

CoY{dBH{Jh,Kh),dBH{Ji,Ki)) = Y{ f f .(s)(-|f - (36) 

i=idRjR 

where J/ = ■ ■ ■, ji,d) G Jh = {jh,i, ■ ■ ■ Jh,d) G ^+, Ki = ■ ■ ■ ,ki,d) ^ and Kh = 

ikh,i ,..., kh,d) G Z([., /i G {1,..., m}. 


Using Eq. (l36l) . Cov{dBH {Jh, Kh), dB^ {Ji, Ki)) can be calculated via two-dimensional wavelet trans¬ 
form of function g{t, s) = —|t — /2 for each z G {1,... , d}- ]E4 h {Jh,K) and E4 h {Ji,K) can 

be calculated similarly. 

When H is known, the covariance matrix of L can be calculated approximately. We denote the 
approximately calculated covariance matrix by G{H), which is a function of H. 

In view of Eq. (l35]) and Eq. (l3^ . G{H) is a positive definite symmetric matrix and differentiable in 
Ef for Ef G (0, l)*^. So G{H)~^ exists and is continuous in EE for EE G (0,1)'^. Now, we can write the 


two-step estimator. 


a = {A'G{Ho)~^A)-^AG{Ho)-^L, 

Hog = a/2 - 1/2, (37) 


where Ho is the OES estimator. 


Theorem 5. The two-step estimator Hog denoted by Eg. (1371) is asymptotically normal, and has the same 
asymptotic covariance as the GLS estimator. 
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Proof: The asymptotic normality follow from Theorem H] 

By Eq. d^ and Eq. (l36l) . it’s easy to check that G{H) is a positive symmetric matrix and differentiahle 

.s p 

in H for H ^ (0,1). And from Theorem |4l for any c G {1,... , d}. Ho —>• H. So we have 

Tc^oo 

G{Ho)-^ ^ G{H)-\ (38) 

Tc^CXD 

Another, G{H) = D'Yar{S)D, where D = log 2 e • diag((Ed^K(Ti, AT)) ^ ,..., {Ed‘^^„{Jm, K)) ^), 
S'= (S(Ji),..., S'(Jm))'. By Eemma|71 Var(S) —)• Tis/nm,c-^o 

Tc^oo ’ 

G{H)^^ j:L/nm,c. (39) 

Tc —^oo 

^Llnm,c is the asymptotic covariance of L. 

Combining (l3^ and (l39l ). the two-step estimator has the same asymptotic covariance as the GES 
estimator. ■ 


IV. Simulation results and discussion 

In this section, we use the estimator Hog to estimate the Hurst parameter of generated fBs. The data 
we test are generated hy the method of circulant embedding of the covariance matrix 1351, l36l. 

a) Selection of parameters: Before the estimation, the octaves J and the number of vanishing 

moments N (or wavelet) must be chosen. Eor octaves J = (ji, • • • ,jd), we choose all the octaves J that 
satisfy Ji < J < Jm, Ji = • • • ,ji,d) is the lower bound of J, Jm = • • • , jm,d) is the lower 

bound of J. Based on prior studies @, Q, for each i G {1,... , d}, jm,,i is chosen the largest possible 
given the data length of dimension i. ji^i is chosen 3 by the minimum mean square error (MSE). It is 
known that there is a bias-variance trade-off for choosing octaves: the selection of small octave increase 
the bias, but decrease the variance of the estimator. The MSE allows the tradeoff between variance and 
bias. 

We choose the classical Daubechies wavelet, which are orthonormal and have compact support. Erom 
Theorem |4l the number of vanishing moments must be chosen N > 2. An increase in the number of 
vanishing moments N comes with an enlargement of the compact support ([— M, M]) Q. In the case 
of finite time, this will lead to the decrease of the number of available wavelet coefficients nj at each 
octaves. So we choose N = 3. 

b) Estimation performance: The estimator Hog for each H is applied to 500 independent copies of 
size 512 X 512 for the two-dimensional case and 100 independent copies of size 256 x 256 x 256 for 
three-dimensional case. All the results are shown in Tabj!] 
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Tabjl] shows that both Hq and Hog have small MSE. For all the estimated H, the means of Ho and 
Hog are both less than H, and the two estimators have the similar bias, but the variance and MSE of 
Hog are less than that of Ho- Besides, we can also find that the smaller the Hi is, the larger the bias and 
MSE are. 

c) Discussion: From our simulation results, we can conclude that both two-step estimator Hog and 
OES estimator Ho perform very well in the sense of small MSE. 

It can be also seen that Hog perform better than Ho because Hog has smaller variance than Ho, which 
is consistent with the theoretical result we have discussed in Section IITI-DI 
Both estimators are not unbiased because 

Elog 2 S{J) / log2E5(J) = log 2 Ed%H{J,K). 

and all the theorems are based on the case of continuous time, but the data we used is discrete, the means 
of both estimators are less than truth parameters as we can see the simulation results. The bias becomes 
larger as the value of H becomes smaller. For Hi > 0.5, the bias of Hi is small and can be ignored. For 
H^ < 0.5, the bias of Hi is obvious but still small. The estimator for Hi > 0.5 performs better than that 
for Hi < 0.5. 


V. Conclusions 

In this paper, we proposed a class of statistical estimators H = {Hi, ..., Hd) for the Hurst parameters 
H = {Hi,..., Hd) of fractional Brownian field via multidimensional wavelet analysis and least squares. 
The proposed estimators are based on the wavelet expansion of fractional Brownian sheet and a regression 
on the log-variance of wavelet coefficient. We also prove that our estimators H are asymptotically normal 
following the approach in 11211 . The convergence rate of H indicates that the estimation accuracy of H 
can be improved by increasing the square root of sample volume y^nf=i The main difficulty of 
our proof is the asymmetric normality of linear combination of S (see Eemma [^, which is achieved 
by the construction of the sequence of random vectors from wavelet coefficients of fBs and 

the function h{X^) that satisfy the condition of Theorem |2] Using the two-step procedure, we realize 
the generalized least squares (GES) estimator, which has the lowest variance of H. We simulated some 
fractional Brownian sheets and used them to validate the two-step estimator. We find that when Hi > 1/2, 
the estimators are efficient, and when Hi < 1/2, there are some bias. 

The computation code of our method is available upon request. 
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TABLE I 

Estimation performance eor Ho and Hog 






Ho 



Hog 


Dimension 


H 

Mean 

Std 

RMSE 

Mean 

Std 

RMSE 

2D 

7Ti 

= 0.3 

0.2733 

0.0297 

0.0399 

0.2706 

0.0144 

0.0327 


H2 

= 0.3 

0.2732 

0.0291 

0.0396 

0.2705 

0.0139 

0.0326 

2D 

Hi 

= 0.5 

0.4897 

0.0301 

0.0318 

0.4911 

0.0150 

0.0174 


H2 

= 0.5 

0.4855 

0.0322 

0.0353 

0.4889 

0.0155 

0.0191 

2D 

Hi 

= 0.8 

0.7915 

0.0323 

0.0334 

0.7958 

0.0165 

0.0170 


H2 

= 0.8 

0.7917 

0.0344 

0.0354 

0.7961 

0.0171 

0.0176 

2D 

Hi 

= 0.3 

0.2731 

0.0303 

0.0405 

0.2708 

0.0140 

0.0324 


H2 

= 0.5 

0.4865 

0.0312 

0.0340 

0.4892 

0.0146 

0.0181 

2D 

Hi 

= 0.3 

0.2736 

0.0324 

0.0418 

0.2719 

0.0153 

0.0320 


H2 

= 0.8 

0.7931 

0.0296 

0.0304 

0.7980 

0.0150 

0.0151 

2D 

7Ti 

= 0.5 

0.4878 

0.0339 

0.0360 

0.4901 

0.0164 

0.0192 


H2 

= 0.8 

0.7911 

0.0337 

0.0348 

0.7971 

0.0152 

0.0155 

2D 

Hi 

= 0.7 

0.6886 

0.0349 

0.0367 

0.6937 

0.0174 

0.0185 


H2 

= 0.8 

0.7902 

0.0334 

0.0347 

0.7959 

0.0169 

0.0174 


Hi 

= 0.6 

0.5929 

0.0171 

0.0184 

0.5943 

0.0085 

0.0102 

3D 

H2 

= 0.7 

0.6942 

0.0172 

0.0181 

0.6973 

0.0080 

0.0084 


Hi 

= 0.8 

0.7955 

0.0152 

0.0158 

0.7970 

0.0087 

0.0091 


Hi 

= 0.7 

0.6941 

0.0186 

0.0195 

0.6967 

0.0095 

0.0100 

3D 

H2 

= 0.8 

0.7931 

0.0181 

0.0193 

0.7959 

0.0089 

0.0098 


Hi 

= 0.9 

0.8962 

0.0167 

0.0171 

0.8980 

0.0105 

0.0106 


The estimation performance lor OLS estimator is given on the left, and tor two-step estimator is given on the right. Std 
denotes standard deviation, RMSE denotes root of MSE. These values are obtained from 500 realizations of generated fBs of 
size 512 X 512, or 100 realizations of generated fBs of size 256 x 256 x 256. It is shown that the means of Ho and Hog are 
both less than H, and the two estimators have the similar bias, but the Std and RMSE of Hog are less than that of Ho (for 

more detail see Section 113. 
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